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Abstract 



Quantum state diffusion (QSD) as a tool to solve quantum-optical master 
\ equations by stochastic simulation can be made several orders of magnitude 

more efficient if states in Hilbert space are represented in a moving basis of 
excited coherent states. The large savings in computer memory and time 
^ ■ are due to the localization property of the QSD equation. We show how the 

O ■ method can be used to compute spectra and give an application to second 



harmonic generation. 
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Stochastic simulation of trajectories of Hilbert-space vectors or quantum trajectories has 
O I proved to be a powerful new method for the solution of master equations in quantum optics. 
^ • The main advantage of stochastic simulation methods over direct numerical solution of the 
■ master equation stems from the fact that far less computer memory is needed to store a 

^ . Hilbert-space vector than to store a density operator. This advantage often far outweighs 
the main disadvantage of stochastic simulation, namely that many trajectories have to be 
added to obtain good statistics. 

There are two main approaches to quantum trajectories in quantum optics: the relative 
state method and quantum state diffusion (QSD) |^. In the relative state method, 
quantum trajectories are conditional on measurement outcomes. Consequently, a single mas- 
ter equation can be "unraveled" into many different stochastic equations, each corresponding 
to a different measurement scheme. For a given problem, one is free to choose the unraveling 
that results in the fastest algorithm. This flexibility is a major strength of the relative state 
method, when the measurement process is the main interaction with the environment. The 
relative state method is limited, however, to problems with few degrees of freedom and to 
relatively small photon numbers. For two or more fleld modes and large photon numbers, 
the number of basis states needed for a numerical representation of a Hilbert-space vector 
can become very large, leading to prohibitive computing times. 

For the QSD method, on the other hand, there is a unique correspondence between master 
equations and stochastic equations @,||. Although it has been shown that, under special 
circumstances, the QSD equation is identical to the relative state equation conditioned on 
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heterodyne detection of photons in general the QSD equation is not derived with respect 
to any specific measurement scheme. The QSD method can be applied with equal ease 
whether or not the main interaction with the environment is a measurement process. For 
an interpretation of the QSD equation, see [Q. 

QSD derives its strength as a computational tool from the localization property of the 
QSD equation p|,[7Hll[]. We have shown [|12| that representing QSD trajectories in a moving 



basis of excited coherent states instead of representing them in the usual number-state basis 
can reduce the number of basis states needed — and thus the computing time — by several 
orders of magnitude. We call the numerical method based on localization quantum state 
diffusion with a moving basis (MQSD). MQSD is particularly advantageous for problems 
with several degrees of freedom and large photon numbers and can be used in regimes 
extending from the highly quantum mechanical (few photons) to the semiclassical (very 
many photons). 

In this paper, we show how MQSD can be used to compute quantum-optical spectra. 
Whereas computing spectra using the relative state method is a well established procedure 
p!3| , quantum state diffusion has not been applied to the practical computation of spectra 
before. Here we show that using MQSD for the computation of spectra leads to the same 
savings as using MQSD for the computation of photon numbers as in |T^, and thus that it 
can be applied where the relative state method would be impractical. 

Many problems in quantum optics can be formulated in terms of a master equation of 
Lindblad form [0 



where p is the system density operator, H is the system Hamiltonian, and the Lm are Lind- 
blad operators which are generally not Hermitian and represent the effect of the environment 
on the system in the Markov approximation. The QSD equation derived from Eq. (|T]) is 
a nonlinear stochastic differential equation for a normalized state vector \ip): 

\d^) = 'h \^lj)dt + Y: {{Ll)L^ - \LlL^ - \{Ll){L^)) mt 

+ J2{Lm-{Lm)) \^)dU. (2) 

m 

The first sum in (0) represents the deterministic drift of the state vector due to the environ- 
ment, and the second sum the random fluctuations. Angular brackets denote the quantum 
expectation (G) = of the operator G in the state lip). The d^m are independent 

complex differential Gaussian random variables satisfying the conditions 

MdU = Md^ndU = , MdCJU = 5nmdt , (3) 

where M denotes the ensemble mean. The density operator is given by the mean over the 
projectors onto the quantum states of the ensemble: 

p = M|7A)(^|. (4) 

If the pure states of the ensemble satisfy the QSD equation (|]), then the density operator 
satisfies the master equation (|lD. 
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A crucial property of the QSD equation is localization. The Schrodinger evolution of an 
isolated system usually produces delocalization or dispersion. According to the Schrodinger 
equation, initially localized wave packets become highly delocalized except under very spe- 
cial circumstances. This effect is also present in the first, Hamiltonian term in the QSD 
equation (^. The environment terms in the QSD equation, however, have the opposite 
effect; they tend to make the wave packet narrower. The theory of localization is treated in 
IQJ^. Numerical simulations show the competition between the delocalizing effect 

of the Hamiltonian and the localizing effect of the environment. The net effect is often a 
very well localized wave packet. 

One important consequence of the localization of quantum trajectories is that, by con- 
tinually changing the basis, it is often possible to reduce the number of basis states needed 
to represent the wave packet by several orders of magnitude. If a wave packet is localized 
about a point {q,p) in phase space far from the origin, it requires many number states \n) 
to represent it. But relatively few excited coherent states \q,p,n) = D{q,p)\n), are needed, 
with corresponding savings in computer storage space and computation time. The operator 
ID{q,p) is the usual coherent state displacement operator [|T^, 



D{q, p) = exp ^ (^pQ - qPj , (5) 
where Q and P are the position and momentum operators. The separation of the represen- 



tation into a classical part {q,p) and a quantum part \q,p,n) is called the moving basis ||T2 
]rU|, the mixed representation. 



or, as m 



The QSD equation can contain both localizing and delocalizing terms. Nonlinear 
terms in the Hamiltonian tend to spread the wave function in phase space, whereas the 
Lindblad terms localize. Accordingly, the width of the wave packets varies along a typical 
trajectory. We use this to reduce the computing time even further by dynamically adjusting 
the number of basis vectors. Details about the implementation of MQSD can be found in 
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We now turn to the computation of optical spectra. The method we use simulates the 



measurement of a spectral component at offset frequency u and is similar in spirit to fTTHlQ 
It involves an auxiliary field mode or filter mode weakly coupled to a system operator c, 
detuned by the frequency u, and weakly damped with damping constant k. The master 
equation describing the system coupled to the filter mode can be written as 

j/ = ^P- ^[Hf, p] + 2K{hpP - \Pbp - \pm) , (6) 

where the system superoperator C is defined as in (|1]), 6 is the annihilation operator for the 
filter mode, and 

Hi = huo Ph + ihe {cP - c^b) (7) 

is the interaction picture Hamiltonian describing the filter mode and its coupling to the 
system with coupling constant e. As in a real experiment, the spectral resolution is limited 
by the filter damping constant n. 
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Assuming that the couphng e is so small that there is a negligible probability for the 
excitation of more than one photon in the filter mode, we can write the total density operator 
as 

p ~ Poo ® |0)(0| + poi ® |0)(1| + pio ® |1)(0| + pn ® |1)(1| , (8) 

where the pij operate in the system Hilbert space, and the operate in the filter Hilbert 
space. By expanding the solution to the master equation in powers of the coupling e, we 
obtain the following approximate equations of motion for the system operators poo and poi 
PPI — notice that these equations do not preserve the trace. 

^poo = Cpoo + 0{e^) , (9) 

^Poi = Cpm + iuj poi + e pooc^ - /tpoi + ^(e^) . (10) 
These equations can be solved formally to yield 

Poo(t) = e^*poo(0)+0(e2) , (11) 

Poi(t) = e f e^(*-*') (poo(t')cO e^^--'')^*-*') dt' + 0{e') . (12) 

J 

We can now obtain an approximation to the Fourier transform of the stationary-state 
time correlation function (c^(0)(i(r))ss for an arbitrary system operator d by computing the 
expectation (b'^dbw) in the stationary regime (by (■ ■ ■)ss we denote the expectation value in 
the stationary state): 

tr [p{0dbb^) = tr [b¥ (poi(t) ® |0)(1|) b^d] = tr [(poi(t) ® \0){0\)d] = tr [dpoi{t)) 

"f/e^(*-*') (poo(t')c^)] e^'--'^)^*-*') dt' (13) 



t 

tr 



{c\t')d{t))e^'''-^^^'-''Ut' , 

where the quantum regression theorem has been used in the last line. In the limit 

t ^ oo, one obtains 



{b^dbb^),,c^e / (c^(0)rf(r))sse(")"rfr . (14) 
Jo 

Similarly, one finds 

POO ^ 

(66^d6)ss^e / (rf(r)c(0))sse(-")^t/r . (15) 
Jo 

By computing the expectation values (b^dbb^)^^ and {bb^db)ss for various choices of the system 
operators c and d using straightforward simulation of the QSD equation, one can thus obtain 
approximations for various spectral densities at the offset frequency ut. For each value of a;, 
a new simulation is needed. 

In the following, we apply the method outlined above to the problem of second harmonic 
generation or frequency doubling. The system consists of two optical modes of frequency Ui 
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and to'2 = which interact in a cavity driven by a coherent external field with frequency 
ijjx and amplitude E. The Hamiltonian in the interaction picture is 



il = ihE{a[ — di) + ih^{a\^d2 — dld^] 



(16) 



where di and 0,2 are the annihilation operators of the two cavity modes, and x describes 
the strength of the nonlinear interaction between them. Damping of the two cavity modes 
is described by the Lindblad operators Li = y/2j^di and L2 = ^27^02. The factors of a/2 
are a consequence of the convention used in the master equation (|^) which differs from the 
usual convention in quantum optics. 

Earlier numerical treat- 



The master equation for this problem first appeared in 22 



ments include [12,23-27]; none of these compute spectra. The main difficulty in a numerical 
treatment of second harmonic generation is due to the large size of the number-state basis 
needed to represent a Hilbert-space vector. If rii number states are needed for mode 1 and n2 
number states are needed for mode 2, then the total number of basis states needed is nin2, 
which can be very large. By using MQSD, and thus exploiting the localization property of 
the QSD equation, this number can be reduced significantly |12[ . We are then in a position 
to compute spectra with moderate numerical effort. 

Of particular interest is the squeezing spectrum in the upconverted mode 0,2 as observed 
in homodyne detection. Squeezing in the upconverted mode in second harmonic generation 
has been experimentally observed by Sizmann et al. [^. The squeezing spectrum Sg^^i^u) 
||ll| is defined in terms of the quadrature operators 



Ae{t) = - d2it)e-'' + dl{t)e 



and 



AAeit) = Ae{t) - {Ae{t)) 



(17) 



(18) 



as 



SeA^) = I672 / dr cos(cur)e-''"(: Aie(0)Aie(r) . 



(19) 



The notation (: ■ ■ ■ :) means normal ordering and time ordering of the operator products, 
and a nonzero value of k results in a limited spectral resolution. A few lines of algebra lead 
to 



Sg^^{uj) = 472 Re 
2k 



dr e 



- 'e*- + e-'^n ie-'''"{d2{T)d2{0)%, + (a5(0)d2(r)). 



(20) 



-2ie 



(«2)ss+ (a2)ss(d2) 



All the quantities in this equation can be easily determined by simulation of the QSD 



equation with the help of Eqs. (14) and ([T5|) . 

In this paper, we do not attempt to reproduce the results of Sizmann et al.; instead, 
we choose a set of parameter values leading to a limit cycle in the corresponding classical 
dynamics. The numerical treatment of this regime has proven to be particularly challenging 
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p3| . Figure |l| shows the resulting spectrum Eq. (|T9D with 6 = 7i/2. The spectral density 
is positive throughout; there is no squeezing in this regime, which agrees with the findings 
in [Q. The stationary expectation values required for the evaluation of Eq. (|20| ) were 
generated by computing time averages in the stationary regime, assuming ergodicity. The 
integration time was T = IOOk^^. Since the filter decay time was chosen to be the 
largest decay time in the system, the statistics resulting from an integration time T = nK~^ 
corresponds roughly to the statistics resulting from adding n trajectories. The computing 
time was about 1 day per data point on a standard Pentium PC running Linux. 

In conclusion, we have shown that the method of quantum state diffusion with a moving 
basis or MQSD can be used to compute quantum optical spectra, with considerable savings 
in computing time and memory compared with earlier methods. 

We thank the EPSRC in the UK and the University of Geneva for essential financial 
support, and G. Alber, N. Gisin, P. Knight, M. Plenio, W. Strunz and P. Zoller for valuable 
communications. 
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FIGURES 

FIG. 1. Homodyne spectrum Eq. ( p!9|) with Q = ■nj'l for the upconverted mode in second 
harmonic generation. The frequency scale on the plot is in units of the decay constant of the 
fundamental mode, 71. The parameters are -E/71 = 20, x/71 = 0-4, 72/71 = !> '^/ti = 0.1, 
e/71 = 10~^, corresponding to a limit cycle regime of the classical dynamics. 
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